/* ------------------------------- Gas-phase ------------------------------- */
Foam::Info << "Create gas-phase mesh for time = " << runTime.timeName() << Foam::endl;

Foam::fvMesh mesh_gas
(
    Foam::IOobject("gas", runTime.timeName(), runTime, Foam::IOobject::MUST_READ)
);

Foam::pointMesh pointMesh_gas(mesh_gas);

// Force calculation of extended edge addressing
{
    const Foam::labelListList& edgeFaces = mesh_gas.edgeFaces();
    const Foam::labelListList& edgeCells = mesh_gas.edgeCells();
    const Foam::labelListList& pointCells = mesh_gas.pointCells();
    Foam::Pout << edgeFaces.size() << " edgeFaces, " << edgeCells.size() << " edgeCells, " << pointCells.size() << " pointCells" << Foam::endl;
}

// Extended stencil of cells (parallel boundaries are automatically handled)
const Foam::extendedCentredCellToCellStencil& addressing = Foam::centredCPCCellToCellStencilObject::New(mesh_gas);

// Initialize stencil storage
Foam::List<Foam::vectorList> sten_pos(mesh_gas.nCells());
Foam::List<Foam::scalarList> sten_ib(mesh_gas.nCells());
Foam::List<Foam::vectorList> sten_val_U(mesh_gas.nCells());
Foam::List<Foam::scalarList> sten_val_p(mesh_gas.nCells());
Foam::List<Foam::scalarList> sten_val_T(mesh_gas.nCells());

// Store coordinates
addressing.collectData(mesh_gas.C(), sten_pos);

/* ------------------------------ Solid-phase ------------------------------ */
Foam::Info << "\nCreate solid-phase mesh for time = " << runTime.timeName() << Foam::endl;

Foam::fvMesh mesh_solid
(
    Foam::IOobject("solid", runTime.timeName(), runTime, Foam::IOobject::MUST_READ)
);

Foam::pointMesh pointMesh_solid(mesh_solid);

// Force calculation of extended edge addressing
{
    const Foam::labelListList& edgeFaces = mesh_solid.edgeFaces();
    const Foam::labelListList& edgeCells = mesh_solid.edgeCells();
    const Foam::labelListList& pointCells = mesh_solid.pointCells();
    Foam::Pout << edgeFaces.size() << " edgeFaces, " << edgeCells.size() << " edgeCells, " << pointCells.size() << " pointCells" << Foam::endl;
}

// Vertex-Cell connectivity (parallel boundaries are considered)
Foam::List<Foam::vectorList> pntAdjCentroid;
collectData_pointCell(pointMesh_solid, mesh_solid.C(), pntAdjCentroid);

/* Overlapping interpolation */
Foam::labelList pnt_cpu;
Foam::labelList pnt_idx;
Foam::vectorList pnt_pos;
Foam::labelList pnt_cIdx;
{
    // Collect gas-phase points over all processors
    Foam::List<Foam::vectorList> pL(Foam::Pstream::nProcs());
    pL[Foam::Pstream::myProcNo()].resize(mesh_gas.points().size());
    for (int i = 0; i < mesh_gas.points().size(); i++)
        pL[Foam::Pstream::myProcNo()][i] = mesh_gas.points()[i];

    Foam::Pstream::gatherList(pL);
    Foam::Pstream::scatterList(pL);

    // Store points within current solid-phase partition
    for (int i = 0; i < pL.size(); i++)
        for (int j = 0; j < pL[i].size(); j++)
        {
            const auto &p = pL[i][j];
            const auto cI = mesh_solid.findCell(p);
            if (cI != -1)
            {
                pnt_cpu.append(i);
                pnt_idx.append(j);
                pnt_pos.append(p);
                pnt_cIdx.append(cI);
            }
        }

    // TODO: check again to cover all corner cases
}

// construct send and recv map
Foam::labelListList sendMap(Foam::Pstream::nProcs());
Foam::labelListList recvMap(Foam::Pstream::nProcs());
Foam::label constructSize = 0;
{
    Foam::labelList nSend(Foam::Pstream::nProcs(), 0);
    for (int i = 0; i < pnt_cpu.size(); i++)
        nSend[pnt_cpu[i]]++;

    for (int i = 0; i < sendMap.size(); i++)
        sendMap[i].setSize(nSend[i]);

    nSend = 0;
    for (int i = 0; i < pnt_cpu.size(); i++)
    {
        const auto procI = pnt_cpu[i];
        sendMap[procI][nSend[procI]++] = i;
    }

    Foam::labelList nRecv(Foam::Pstream::nProcs(), 0);
    Foam::Pstream::exchangeSizes(sendMap, nRecv);

    for (int i = 0; i < recvMap.size(); i++)
        recvMap[i].setSize(nRecv[i]);

    // Construct with my own elements first
    for(int i = 0; i < recvMap[Foam::Pstream::myProcNo()].size(); i++)
        recvMap[Foam::Pstream::myProcNo()][i] = constructSize++;

    // Construct from other processors
    for(int i = 0; i < recvMap.size(); i++)
    {
        if (i != Foam::Pstream::myProcNo())
        {
            for(int j = 0; j < recvMap[i].size(); j++)
                recvMap[i][j] = constructSize++;
        }
    }
}

Foam::mapDistribute map(constructSize, Foam::move(sendMap), Foam::move(recvMap));
